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1.  OVERVIEW 

This  Final  Technical  Report  constitutes  the  completion  of  the  PETTT  Pre-Planned  Project  PP- 
SAS-KY06-001-P3. 

The  development  of  predictive  capabilities  for  modeling  spacecraft  environments  in  the  near¬ 
space  regime  for  re-entry  applications  and  for  satellite  trajectories  presents  significant  challenges 
due  to  the  multi-physics,  multi-scale  nature  of  the  flows  of  interest.  A  challenge  directly  related 
to  the  multi-scale  flows  over  space  vehicles  is  accurate  numerical  prediction  of  flows  that  contain 
most  or  all  gas  dynamic  regimes  from  continuum  to  rarefied  to  free  molecular.  To  address  this 
problem,  the  goal  of  this  project  is  the  improvement  of  gas  dynamics  solvers  that  are  currently 
used  across  several  DoD  applications.  Significant  steps  were  taken  toward  the  goal  by  achieving 
the  following  objectives 

(1)  The  implementation  of  a  new  BGK-type  model  with  velocity-dependent  collision 
frequency  and  an  increased  fidelity  in  relaxation  of  moments, 

(2)  The  development  of  methods  for  fast  evaluation  of  the  Boltzmann  collision  integral 
based  on  the  stochastic  integration, 

(3)  The  validation  of  the  improved  physical  accuracy  of  the  new  modules. 

1.1  Summary  of  the  Accomplishments 

1.1.1  Objective  1.  The  implementation  of  a  new  BGK-type  model  with  velocity-dependent 
collision  frequency  and  an  increased  fidelity  in  relaxation  of  moments. 

The  following  major  tasks  were  completed  to  achieve  the  first  objective  of  the  proposal. 

(1)  A  stand-alone  library,  DGVlib,  was  developed  to  encompass  the  following  software 
tools: 

a.  Subroutines  for  nodal  discontinuous  Galerkin  (DG)  discretizations  of  the  kinetic 
equations  in  the  velocity  variables  on  uniform  and  octree  meshes; 

b.  Subroutines  for  evaluating  the  collision  operator  using  BGK-type  model  with 
velocity  dependent  collision  frequency; 

c.  Subroutines  for  evaluating  Boltzmann  collision  operator  using  nodal  DG 
discretizations  on  uniform  grids  and  subroutines  to  pre-compute  values  of  the 
collision  kernel. 

(2)  The  BGK  model  with  velocity  dependent  collision  frequency  was  generalized  to  allow 
for  a  large  number  of  moments  to  be  enforced  approximately.  Subroutines  implementing 
the  extended  model  were  added  to  the  library. 

(3)  Techniques  for  automatic  selection  of  the  kinetic  collision  model  were  developed  and 
implemented  in  the  library.  The  evaluation  of  the  kinetic  collision  operator  can  be 
accomplished  by  a  call  of  a  single  subroutine  which  will  choose  an  appropriate  kinetic 
model  based  on  the  norm  of  deviation  of  the  solution  from  the  continuum  state. 

(4)  The  DGVlib  library  has  been  merged  with  a  number  of  drivers  implemented  in  Fortran. 
Specifically,  it  was  merged  with  a  high  order  multi-step  time  integration  driver  in  zero 
spatial  dimensions  and  with  a  high  order  Runge-Kutta  DG  driver  in  one  spatial 
dimension  that  were  available  to  the  PI.  The  library  was  merged  with  a  second  order 
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accurate  finite  volume  two  dimensional  solver  SMOKE  available  to  the  co-PI 
Gimelshein.  The  software  that  resulted  from  the  merger  was  used  to  perform  validation 
of  the  library  accuracy  described  in  Objective  3.  The  library  was  also  merged  with  a 
C++  driver,  the  Unified  Flow  Solver. 

(5)  A  detailed  user  guide  was  developed  to  help  the  user  utilize  DGVlib. 

1.1.2  Objective  2.  The  development  of  methods  for  fast  evaluation  of  the  Boltzmann 
collision  integral  based  on  the  stochastic  integration. 

For  the  second  objective,  a  new  numerical  method  was  designed  for  evaluating  nodal-DG 
discretizations  of  the  Boltzmann  collision  operator  on  uniform  grids  based  on  Korobov 
integration.  The  method  requires  pre-computing  and  storing  values  of  the  kernel  of  the  collision 
operator  at  quasi- stochastic  nodes.  Subroutines  for  computing  and  storing  values  of  the  kernel 
and  subroutines  for  evaluating  the  collision  integral  using  quasi- stochastic  Korobov  quadratures 
were  developed  and  added  to  the  library.  Accuracy  tests  of  the  new  method  were  performed. 

1.1.3  Objective  3.  Validation  of  the  improved  physical  accuracy  of  the  new  modules. 

Solutions  to  the  following  problems  were  computed  using  the  new  models  of  DGVlib  library. 

(1)  Problems  of  spatially  homogeneous  relaxation  were  solved  to  study  in  detail  the 
properties  of  the  new  models  of  collision  integral. 

(2)  One  dimensional  problems  of  normal  shock  wave  were  solved  for  Mach  numbers  1.55, 
3.00,  3.8,  6.5,  and  9.0. 

(3)  Two  dimensional  super-sonic  flow  past  cylinder  was  computed. 

Zero-dimensional  solutions  were  compared  to  solutions  to  the  full  Boltzmann  equation.  One  and 
two  dimensional  solutions  were  compared  to  DSMC  solutions  computed  by  an  established  code 

SMILE.  Solutions  computed  by  new  models  showed  a  significant  improvement  in 
predicting  temperature  profiles  of  the  flows  and  improvement  in  predicting  density  profiles 
of  the  flows  as  compared  to  the  classical  ES-BGK  model. 

1.2  Areas  of  Applicability  of  the  Proposed  Software 

The  developed  software  capability  is  expected  to  provide  essential  elements  for  the  integration 
and  generalization  of  kinetic  and  hybrid  solvers,  both  currently  used  and  to  be  developed  in 
DoD.  It  primarily  targets  various  AF  applications  where  flow  non-equilibrium  is  important,  such 
as  hypersonic  flight  at  moderate  to  high  altitudes  and  rocket  propulsion  with  plume  -  atmosphere 
interaction,  both  for  conventional  and  miniaturized  thrusters,  as  well  as  transient  rarefied  flows. 
DGVlib,  made  available  to  DoD  researchers  and  engineers,  will  allow  for  expanding  the  area  of 
applicability  of  the  existing  fluid  dynamics  codes,  and  stimulate  the  development  of  state-of-the- 
art  hybrid  codes  applicable  to  all  flow  regimes  from  continuum  to  free  molecular. 

1.3  Potential  Impact  of  the  Project 

The  Air  Force  is  focused  on  the  development  of  advanced  space  technologies  for  more  effective, 
more  affordable  warfighter  missions.  To  make  progress  in  this  area,  it  is  necessary  to  develop 
predictive  capabilities  for  modeling  spacecraft  environments  in  the  near-space  regime  for 
reentry 
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applications  as  well  as  for  satellite  trajectories.  Other  DoD  programs  that  rely  critically  on 
technologies  in  this  flight  regime  are  AF's  global  strike  weapons  and  boost  glide  vehicles,  as  well 
as  MDA's  exo atmospheric  anti-ballistic  missiles  programs.  The  availability  of  efficient  and 
physically  accurate  predictive  tools  for  gas  flow  modeling  in  a  wide  range  of  flow  regimes  is 
critical  for  the  evaluation  of  local  and  overall  aerothermodynamic  loads  and  the  detailed  analysis 
of  existing  and  future  propulsion  systems.  Computer  model  based  tools  are  also  one  of  the 
central  components  of  the  space  situational  awareness,  where  the  comprehensive  knowledge  of 
space  objects  and  the  ability  to  track,  understand  and  predict  their  future  location  becomes  the 
key. 

The  developed  software  library  addressed  several  improvements  in  the  existing  solvers. 

•  The  implemented  BGK  model  collision  operator  with  velocity  dependent  collision 
frequency  was  shown  to  increase  physical  accuracy  in  predicting  temperature  profiles  of 
non-continuum  flows.  With  the  availability  of  this  module  a  user  gains  additional 
flexibility  on  what  properties  of  the  solution  may  be  enforced  in  the  kinetic  models  at 
moderate  additional  computational  costs. 

•  The  implemented  module  performing  integration  of  the  collision  integral  using  Korobov 
nets  computes  the  Boltzmann  collision  operator  faster  than  full  deterministic  integral 
evaluation. 

•  The  modules  that  implement  fully  deterministic  solution  of  the  Boltzmann  equation  can 
be  conveniently  used  to  provide  accurate  solutions,  although  these  solutions  will  not  be 
efficient. 

•  The  implemented  algorithms  for  choosing  a  different  kinetic  model  based  on  the  norm 
of  deviation  of  the  solution  from  continuum  provide  a  convenient  tool  to  use  in  kinetic 
simulations  when  the  most  efficient  model  can  be  used  to  compute  the  collision  operator 
at  a  point. 

Overall,  the  implemented  software  has  high  potential  to  improve  gas  dynamics  solvers  that  are 
currently  used  across  several  DOD  platforms.  Through  that,  the  implemented  software  has 
potential  to  address  the  following  problems  pertaining  to  the  Space  and  Astrophysical  Sciences 
area:  (1)  modeling  of  near-space  vehicle  aerothermodynamic s  for  flight  conditions  where  there 
are  both  continuum  and  rarefied  flow  regions,  (2)  spacecraft  thrusters  that  exhaust  high  density 
plumes  into  rarefied  atmosphere,  (3)  gas  flows  in  microscale  onboard  devices.  The  completed 
project  benefits  the  following  DOD  research  topics:  MDA1 2-022  “Miniature  extendable  nozzles 
or  actuating  nozzles  for  improved  ISP  of  DACS  thrusters”;  AF 14 1-089  “Electric  propulsion  of 
orbit  transfer”,  and  MDA 12-009  “Fast-running  physics-based  models  for  intercept  debris  aero- 
heating  and  aero -thermal  demise”. 

2.  ACCOMPLISHMENTS 

A  software  library,  DGVlib,  has  been  developed  to  encompass  a  set  of  tools  that  can  be  used  to 
solve  kinetic  equations  of  gas  dynamics.  The  name  of  the  library,  DGVlib,  stands  for  nodal- 
Discontinuous  Galerkin  discretizations  in  the  Velocity  variable  (nDGV).  The  library  includes 
subroutines  to  perform  nDGV  discretizations  in  the  three  dimensional  velocity  space  including 
capabilities  for  non-uniform  octree  discretizations,  subroutines  for  evaluating  collision  operators 
in  gas-kinetic  models,  including  the  BGK,  ES-BGK  and  the  new  velocity-dependent  BGK 
model 
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with  velocity  dependent  collision  frequency.  The  library  includes  subroutines  to  evaluate  the 
Boltzmann  collision  integral  using  Korobov  quadratures.  DGVlib  library  was  merged  with  a  high 
order  multiple  time  stepping  driver  in  zero  spatially  dimensions,  a  high  order  Runge-Kutta 
discontinuous  Galerkin  driver  in  one  spatial  dimension,  and  with  the  second  order  accurate  finite 
volume  solver  SMOKE  in  two  spatial  dimensions.  The  resulting  codes  were  used  to  perform 
benchmark  simulations  which  include  problems  of  spatially  homogeneous  relaxation,  normal 
shock  wave  and  super-sonic  flow  past  cylinder. 

2.1  Key  Components  of  DGVlib 

The  new  software  library  has  the  following  key  components: 

1 .  A  collection  of  subroutines  to  build  elements,  nodes,  and  basis  functions  of  nodal 
discontinuous  Galerkin  discretizations  of  functions  in  three  dimensional  velocity  space. 

2.  Subroutines  for  evaluation  of  the  collision  operators  in  gas  kinetic  models,  including  the 
BGK,  the  ES-BGK,  and  the  new  velocity-dependent  BGK  model  with  velocity 
dependent  collision  frequency. 

3.  Subroutines  for  the  evaluation  of  the  Boltzmann  collision  operator  using  Korobov 
quadratures. 

4.  Interfaces  to  call  collision  operator  subroutines  from  other  solvers  in  zero,  one  and  two 
spatial  dimensions. 

2.2  DGVlib  Components  and  Functionality. 

The  overall  structure  of  the  DGVlib  is  shown  in  Figure  1.  Key  capabilities  and  main  features  of 
the  library  are  described  below. 

2.2.1  Subroutines  Implementing  Nodal  Discontinuous  Galerkin  Discretizations  in  the  Velocity 
Variable. 

DGVlib  includes  subroutines  to  create  rectangular  meshes  in  3D  velocity  space.  It  also  includes 
subroutines  to  create  octree  non-uniform  meshes  in  3D  velocity  space;  however,  octree  meshes 
are  not  supported  by  the  subroutines  evaluating  the  collision  operators.  At  the  core  of  the  library 
are  subroutines  implementing  high  order  nodal  discontinuous  Galerkin  discretization  in  the 
velocity  space  (nDGV)  [JJ.  Subroutines  for  interpolating  nDGV  solutions  from  the  Galerkin 
coefficients  are  also  included.  The  library  maintains  two  nDGV  discretizations,  referred  as  the 
primary  and  the  secondary  meshes  in  Figure  1.  nDGV  discretizations  are  functionally  very 
similar  to,  albeit  more  powerful  than,  discrete  ordinate  discretizations  of  kinetic  equations  which 
are  commonly  used  for  simulations  of  rarefied  gas  dynamics.  As  a  result,  DGVlib  can  be  merged 
relatively  easily  with  kinetic  solvers  that  are  based  on  the  discrete  ordinate  approach.  A  folder 
view  of  DGVlib  merged  with  a  zero  dimensional  solver  is  shown  in  Figure  2. 

2.2.2  Subroutines  Implementing  Collision  Operators  in  Gas-Kinetic  Equations. 

DGVlib  includes  several  subroutines  that  implement  evaluation  of  the  collision  operators  in  gas 
kinetic  equations.  The  following  models  are  supported:  the  BGK  model,  the  ES-BGK  model,  and 
the  new  BGK  model  with  the  velocity  dependent  collision  frequency  [2].  In  addition,  the  library 
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implements  evaluation  of  the  full  Boltzmann  collision  operator  using  the  formulation  of  [I] 
directly  and  using  Korobov  quadratures.  The  BGK  model  with  velocity-dependent  collision 
frequency  relies  on  evaluation  of  the  full  collision  operator  to  determine  the  correct  relaxation 
rates  in  the  solution.  All  models  are  supported  on  the  primary  mesh.  Also,  deterministic 
evaluation  of  the  full  collision  operator  is  supported  on  the  secondary  mesh.  Model  collision 
operators  are  not  yet  supported  on  the  secondary  mesh,  since  there  has  been  no  user  demand  for 
this  feature.  A  conglomerate  subroutine  is  provided  for  the  evaluation  of  the  collision  operator  on 
the  primary  mesh,  called  UniversalCollisionOperatorDGV  that  automatically  chooses 
the  model  based  on  the  degree  of  deviation  of  the  solution  from  the  local  Maxwellian 
distribution. 


Figure  1.  Functionality  of  DGVIib  within  a  kinetic  solver. 


2.2.3  The  Use  of  DGVIib  with  OD,  ID,  and  2D  Solvers. 

The  library  is  to  be  called  from  an  external  driver.  First,  DGVIib  subroutines  for  creating  nDGV 
discretizations  are  called  to  construct  the  discretization  in  the  velocity  variable.  Then,  the  arrays 
that  implement  the  velocity  discretization  are  passed  to  the  driver  to  construct  the  full  spatial  and 
velocity  discretization.  In  the  case  of  ID  and  2D  applications,  discretizations  in  the  spatial 
variables  are  provided  by  the  driver.  Currently,  the  complete  spatial  and  velocity  discretization  is 
obtained  by  taking  the  full  tensor  product.  As  a  result,  all  spatial  cells  have  the  same  velocity 
discretization.  To  invoke  kinetic  collision  operators  in  ID  and  2D  applications,  subroutines  are 
provided  that  evaluate  the  collision  operator  over  a  collection  of  spatial  cells  referenced  by  a 
single  index.  The  library  maintains  data  structures  to  store  auxiliary  information  about  the 
solution  at  each  spatial  cell.  In  particular,  the  library  stores  relaxation  rates  computed  for  each 
spatial  cell  that  are  used  in  the  BGK  model  with  velocity  dependent  collision  frequency.  It  also 
keeps  track  of  when  the  rates  have  to  be  updated. 
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Figure  2.  A  folder  view  of  DGVIib  library  merged  with  a  zero  dimensional  driver. 

2.3  Mathematical  Models  Utilized  in  DGVIib 

DGVIib  library  encompasses  subroutines  implementing  the  recently  proposed  BGK-type 
(Bhatnagar-Gross-Krook)  model  with  velocity-dependent  collision  frequency  and  increased 
fidelity  in  relaxation  of  moments  and  subroutines  that  implement  integration  of  the  Boltzmann 
collision  integral  using  Korobov  nodes.  In  addition,  subroutines  for  evaluating  the  classical  BGK 
model  and  the  ellipsoidal-statistical  BGK  model  and  subroutines  for  evaluating  the  full 
Boltzmann  collision  operator  using  the  formulation  of  [I]  are  also  included  in  the  library.  The 
following  sections  briefly  summarize  the  mathematical  models  that  are  employed  in  the  first  two 
approaches.  Additional  discussion  and  also  detail  on  the  other  implemented  models  can  be  found 
in  Appendix  A  to  this  report. 

2.3.1  Dimensionless  Reduction  of  the  Velocity  Variable 

Gas  dynamic  constants  vary  greatly  in  scale.  Caution  should  be  exercised  when  these  quantities 
are  combined  in  a  numerical  method  in  order  to  avoid  accumulation  of  the  round-off  errors. 
However,  because  of  the  large  number  of  integration  points,  accumulation  of  round-off  errors  is 
difficult  to  avoid  in  evaluation  of  the  collision  operator.  The  use  of  dimensionless  reduction 
when  all  the  key  quantities  are  of  the  order  of  unit  can  help  avoid  many  obstacles  caused  by  the 
round-off  errors.  DGVIib  uses  the  dimensionless  reduction  that  can  be  found  in  Chapter  3  in  [3]. 
Below,  the  key  relations  are  very  briefly  summarized. 
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Let  x,  v  be  the  conventional  dimensional  variables  of  time,  space  and  velocity.  Their 
dimensionless  counterparts  are  defined  by 


t 

1  ~~  T’ 


x 

1: 


V 


V  = 


Cn 


or  t  =  £T,  x  =  xL ,  v  =  vCc 


Where  ^oc  is  some  reference  temperature  and  L  is  some  reference  length  scale.  The  reference 
velocity  Coo  is  conventionally  defined  by  Coo  —  V^RT^  wherc  R  is  the  gas  constant.  The 
reference  time  scale  T  is  defined  by 


T  =  L/C0 o 


The  dimensionless  density,  bulk  velocity  and  temperature  are  defined  by 


n(t,  x ) 


f(t,x,v)dv 


n(t,  x)u(t,  x ) 


vf(t,  x,  v)dv 


n(t,x)T(t,x)  :=  -  [  (v  —  u)2  f(t,  x,v)dv 

3  JR  3 

Here  the  dimensionless  velocity  distribution  function  f  x>  v )  is  given  by 

T  3/^3  ^ 

f(t,x,v), 

where  N  is  the  total  number  of  molecules  in  the  volume  of  gas  L  ' .  Mathematical  models 
presented  in  this  section  are  given  in  the  dimensionless  variables. 

2.3.2  BGK  Model  with  Velocity-Dependent  Collision  Frequency 

This  section  briefly  describes  the  BGK-type  models  with  the  velocity  dependent  collision 
frequency  that  are  implemented  in  the  subroutine  EvalColVelES.  Two  models  are  used  in  this 
approach.  The  first  model  will  attempt  to  enforce  exactly  the  relaxation  rates  for  the  selected 
group  of  moments.  In  this  approach  the  number  of  enforced  moments  and  the  number  of  basis 
functions  in  the  representation  of  the  collision  frequency  have  to  be  the  same.  A  linear  system  of 
equations  is  solved  at  every  time  step  to  determine  the  coefficients  in  the  representation  of  the 
collision  frequency.  In  the  second  model,  the  number  of  enforced  moments  is  greater  than  the 
number  of  basis  functions.  A  linear  least  squares  problem  is  solved  at  every  time  step  to 
determine  the  coefficients  of  the  collision  frequency.  The  two  models  are  summarized  below. 

2.3.3  BGK-Type  Model  with  Exact  Enforcement  of  Moments 

Let  llJt  (' u )  ’  ~  0, ...  ,k  be  a  collection  of  linearly  independent  functions  in  velocity  space.  To 
allow  the  BGK  model  freedom  to  accommodate  the  specified  relaxation  frequencies  for  selected 
solution  moments  one  assumes  that  the  collision  frequency  has  the  form 
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k 

v{t}u)  =  y ^Ci(t)if>i(u). 

i= 0 

In  general,  it  is  best  to  select  functions  ^i\u)  so  that  the  projection  mapping  between  the  kernels 
of  the  enforced  moments  and  the  functions  is  invertible  and  is  well-conditioned  in  a 

weighted  L2-norm  with  f M{t ■,  u)  serving  as  the  weight.  In  simulations  presented  in  the  Section 
3  the  following  functions  were  used: 

0o(^)  =  1,  01  (u)  =  Ui,  02(u)  =  U2,  03(u)  =  u3.  04 ( u)  =  u\+u\+U%, 

00m)  =  Pi-3(u i),  i  =  5, . . . ,  8, 

where  ^  Oi)  are  the  Legendre  polynomials  of  degree  l.  It  is  assumed  that  the  coefficient 

standing  with  ^’o(w)  has  the  form  co  +  z/hgk.  Thus  one  seeks  an  addition  to  the  collision 
frequency  of  the  classical  BGK  model  that  enforces  relaxation  speeds  of  selected  moments.  With 
this  assumption,  the  new  model  may  be  considered  as  a  generalization  of  the  classical  BGK 
model. 

The  coefficients  Ci  (^)>  *  —  1:  .  .  .  ,  k  arc  determined  from  the  condition  that  moments  fM  of 
the  solution  relax  with  the  prescribed  relaxation  frequencies  Uip.  Specifically,  the  following 
system  is  solved  at  each  time  step: 

k  „ 

Mt)  -  Z'BGK ){f™  ~U)=y2  Ci(t)  /  00m)(/m(0  u)  -  f(t,  u))<p(u)  du  . 

i= o  JR3 

The  relaxation  frequencies  are  determined  by  evaluating  the  full  Boltzmann  collision  integral 
using  the  formula: 


G00 


/ 


v> 


M*)  -  f¥(t) 


Where  is  the  Galerkin  projection  of  the  Boltzmann  collision  operator  up  to  coefficients 

arising  in  dimensionless  reduction.  The  kernels  <p(v)  run  over  the  set  •  •  •  5  [ fkip )  that 

determine  which  properties  are  enforced  in  the  solution. 

2.3.4  Extended  Model  for  Enforcing  Relaxation  Rates 

The  model  with  exact  enforcement  of  moments  has  one  undesired  property,  namely  that  the 
kinetic  models  can  become  stiff  (more  stiff  than  the  classical  models)  if  relaxation  rates  are 
enforced  for  higher  order  moments.  To  reduce  the  stiffness  one  can  change  the  way  in  which  the 
relaxation  rates  of  moments  are  enforced. 
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In  the  model  of  the  previous  section,  the  number  of  functions  ^kip)  and  the  number  of  enforced 

moments  U(t)  is  the  same,  resulting  in  a  linear  system  of  equations  with  a  square  matrix.  In 
principle,  one  may  desire  to  enforce  relaxation  rates  for  a  larger  number  of  moments  than  the 

number  of  function  r(Pi(u).  in  this  case,  the  linear  system  of  the  previous  section  is 
overdetermined,  that  is,  has  a  rectangular  matrix. 


It  will  be  assumed  that  the  collision  frequency  v  is  expressed  in 
and  corresponding  coefficients  cfc(A  x),  k  =  1 ..  .  Nf^  as 


terms  of  basis  functions  ^k{v) 


Nh 

u(t,x,v)  =  yjCk(t)Mv)- 

k= 1 


For  convenience  of  later  use,  a  vector  of  parameters,  c  ,  is  defined  as 

c  =  (cq.  Co ;  •  •  • ;  C/V(, ) 

For  any  function  ^(y),  the  corresponding  (generalized)  moment  can  be  defined  as 


U(t,  x) 


(f>(v)f(t,  x,  v)  dv  . 


Further,  given  any  selection  of  moments  (Piyv)  for )  —  f  •  •  •  Nc,  one  can  obtain  the  rate  of 
evolution  (due  to  collisions)  of  any  moment  based  on  a  weak  form  of  the  Boltzmann  equation  as, 


collision 


where  the  right  hand  side  of  the  above  equation  can  be  conveniently  evaluated  using  certain 
symmetry  properties  of  the  collision  operator.  These  properties  ensure  that  the  above  right  hand 
side  term  is  zero  for  collision  invariants. 


Similarly,  the  rate  of  evolution  obtained  from  the  model  approximations  of  the  collision  operator 
can  be  expressed  as, 

=  f  x,  v)(f0(t,  x,  v)  -  f0(t,  x,  v))  dv. 

Jus 

Based  on  the  assumed  functional  form  of  v,  in  terms  of  basis  functions  i>k{v)  ,  the  rate  of 
evolution  of  the  generalized  moment  can  be  expanded  as, 
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where 


Nb 

,  k=l 


EG- 


1 


fo(t,x,v))dv 


G 


ik 


(pi(v)ipk(v)(fo(t:  x.  v)  -  fo(t,  X,  v ))  dv  . 


In  order  to  determine  the  unknown  coefficients  °k v ' '  x),  one  defines  an  error  measure  e  and 
seeks  those  coefficients  that  minimize  the  error.  The  error  function  e  is  chosen  to  be  in  the 
following  form: 


e 


i—1 


i—  1 


where  ^ c  denotes  the  number  of  constraints  on  evolution  of  moments  that  we  choose  to 
consider.  One  can  see  that  the  minimization  of  e  is  a  classical  linear  least  squares  problem  to 
which  standard  approaches  can  be  applied.  In  particular,  if  the  matrix  G  is  full  rank  then  the 
solution  to  the  problem  is  given  by 


c  =  [GTG\~1GrMlt. 

If  a  matrix  G  is  ill-conditioned  or  even  not  full  rank,  the  problem  can  be  regularized  using  the 
singular  value  decomposition. 

2.3.5  Implementation  of  Korobov  Integration  for  Evaluating  the  Collision  Operator 

Stochastic  evaluation  of  the  Boltzmann  collision  integral  has  been  the  most  successful  approach 
for  obtaining  solutions  in  multidimensional  applications  [4]  Stochastic  methods  are,  indeed,  the 
only  well-developed  methods  that  beat  the  "curse  of  dimensionality"  in  multidimensional 
integration.  A  family  of  efficient  quadratures  is  due  to  Korobov  [5],  who  proposed  the  use  of 
decimal  parts  of  fractions  as  Monte-Carlo  integration  points  and  showed  that  the  resulting  quasi¬ 
stochastic  multidimensional  quadratures  converge  fast  in  the  case  of  smooth  periodic  functions. 

A  description  of  the  Korobov  quadratures  implemented  in  DGVlib  can  be  found  in  Chapter  3 
Section  24  of  [5].  This  section  briefly  summarizes  the  algorithm  implemented  in  DGVlib. 

Let  a  multidimensional  function  h(xi  j  •  ■  •  ;  xs)  be  periodic  with  period  one  on  a  5-dimensional 
cube  ^  ~  ([O’  1])  .  Let  -  -  •  >  xs)  have  an  absolutely  converging  Fourier  series  on  D  with 

Fourier  coefficients  ,  .  .  .  ,  ms)  satisfying  the  inequality 
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C(m i, . .  . ,  ms)  < 


C 


(mi  ■  ■  •  ms)c 


,  where  m  =  max(l,  |m|),  for  some  a  >  1. 


Consider  the  quadrature  formula 


f  •••  f  h{x1,...,xs)dx1...dxa  =  -  — . . . ,  {  —  -  Rp[h\, 

J  0  J  0  P  P  P 


where  p  is  a  large  number,  1  ~  •  •  •  ?  s,  are  some  whole  numbers  relatively  prime  with  p, 

and  {  '  }  stands  for  the  fractional  (decimal)  part  of  the  number.  The  following  estimate  holds  for 
the  truncation  error: 


Rp\h]  |  <  Ci(a,  s ) 


In  asp 

pa 


It  follows  that  the  convergence  rate  a  is  proportional  to  the  smoothness  of  the  function. 
Specifically,  a  higher  smoothness  of  the  function  h  implies  the  faster  convergence  rate  of  the 
quadrature  rule. 


Subroutines  of  DGVlib  implement  the  Korobov  integration  of  the  Galerkin  projection  of  the 
Boltzmann  collision  integral  (see  Appendix  A  and  also  [1]) 


N2d2Cr 


H3p 


LG 


f(t,  x,  v)f(t ,  x,  u)A(v ,  u ;  (f)J)dudv 


N2d2Coo 

LG 


p 

EE  f{t,  X,  vk)f{t ,  x,  uk)A(vk,  Uk]  ftp), 

k= l 


where 


(AV)2 

A{vk,uk\(j? )  = - A{vk,uki(jA) 

p 

and  V  —  [uLi  ur\  x  \vl,  l'R  x  [wl,  wr]  [s  the  velocity  domain,  and  AV  is  the  volume  of 

— >  — t 

the  domain  V.  The  Korobov  nodes  vki  11  k  are  obtained  by  scaling,  i.e., 

-  ( aik  a2k  a3k  \ 

Vk  =  - [UR  ~  UL)  +  UL,  - [Vr  ~  VL)  +  VL,  - [WR  -  WL )  +  WL 

V  p  v  p  ) 

( d’4: k  ,  \  CL§k  .  \  CLftk  ,  \  \ 

uk  =  - [UR  -  UL)  +  UL,  - [VR  ~  VL)  +  VL.  - [Wr  ~  WL)  +  WL 

V  p  v  p  ) 
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The  summation  kernel  ^  (  l'k  ’  Ukl  is  pre-computed  and  stored  for  future  use. 

2.4  Numerical  Accuracy  and  Validation 

2.4.1  Simulations  in  Zero  Spatial  Dimensions  Using  the  Improved  Velocity-Dependent 
Collision  Frequency 

The  Pis  developed  an  improved  model  for  enforcing  relaxation  rates  using  a  velocity-dependent 
collision  frequency.  The  mathematical  approach  to  the  new  velocity-dependent  collision 
frequency  is  described  in  [2,6].  The  new  approach  has  been  implemented  in  the  DGVlib  library. 
This  section  summarizes  the  advantages  of  the  new  model  in  the  solution  to  the  problem  of 
spatially  homogeneous  relaxation  of  two  Maxwellian  streams. 

The  two  Maxwellian  streams  in  the  considered  solutions  have  dimensionless  densities,  bulk 
velocities  and  temperatures,  rij  =  1.609,  ui=  (1.228;  0;  0),  7)  =  0:043  and  n2  =  6.011,  u2  =  (0.329; 
0;  0),  Ti  =  0.603.  The  reference  number  density  is  Nref  =  10E+20  and  the  reference  temperature 
is  Tref=  7000  K.  These  values  of  macroparameters  correspond  to  downstream  and  upstream 
conditions  of  a  normal  shock  wave  with  Mach  number  6.5. 

In  Figure  3,  solutions  obtained  using  the  BGK  model  with  the  improved  velocity-dependent 
collision  frequency,  the  BGK  model  with  original- velocity  dependent  collision  frequency  and  the 
classical  BGK,  ES-BGK  and  the  Shakhov  models  are  compared.  The  difference  between  the 
original  and  the  improved  velocity  dependent  collision  frequencies  is  due  to  the  way  the 
relaxation  rates  for  moments  are  enforced.  In  the  original  model,  the  number  of  enforced 
relaxation  rates  is  equal  to  the  number  of  the  basis  functions  in  the  representation  of  the  velocity 
dependent  collision  frequency.  In  particular,  to  enforce  relaxation  rates  of  high  order  moments, 
higher  order  polynomial  basis  functions  are  used.  The  use  of  high  order  polynomials  results  in 
strong  perturbations  of  moments  that  are  not  enforced.  An  example  of  such  perturbation  can  be 
seen  in  Figure  3(h).  In  this  figure,  relaxation  of  moments  with  kernels  (ui-vj ,  i=l,2,  are  shown. 
Here  v„  i=l,2,3  are  components  of  gas  bulk  velocity.  These  moments  were  not  enforced  in  the 
simulations  and  appear  to  be  strongly  perturbed  in  the  original  model  when  high  order  moments 
are  enforced.  The  improved  model  allows  the  number  of  enforced  moments  to  be  larger  than  the 
number  of  basis  functions.  In  the  new  model,  the  coefficients  of  the  collision  frequency  are 
determined  by  solving  a  linear  least  squares  problem  using  singular  value  decomposition  and 
regularization.  In  Figure  3(g),  relaxation  of  moments  (uj-Vif  \  i=l,2,  is  shown  in  solutions 
computed  by  the  improved  model.  The  case  marked  by  m=5  corresponds  to  solutions  in  which 
moments  with  kernels  (ui-vi)m,  i=l,2,  m=3,4,5  are  enforced  while  only  polynomials  up  to  second 
degree  were  used  in  the  representation  of  the  collision  frequency.  It  can  be  seen  that  the  new 
model  results  in  much  smaller  perturbations  of  moments  that  are  not  enforced  as  compared  to  the 
original  model.  In  Figure  3(d — f)  relaxation  of  the  moments  (urv,)4,  i=l,2,  is  presented  for 
solutions  computed  by  the  new  model  (d),  the  original  model  (e),  and  the  classical  models  (f).  In 
the  original  model  only  moments  (u,-v,/H,  i=l,  m=2,3,4,5  are  enforced  and  high  order 
polynomials  are  used  in  the  representation  of  the  collision  frequency  to  match  the  number  of 
enforced  moments.  The  results  suggest  that  the  new  model  approximates  relaxation  of  high  order 
moments  comparable  to  the  original  model  and  better  than  the  classical  models.  Nevertheless,  as 
the  number  of  enforced  moments  increases  in  the  new  model  and  the  number  of  the  basis 
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functions  is  kept  the  same,  it  becomes  difficult  to  achieve  accurate  rates  for  the  moments  because 
the  residuals  in  the  linear  least  squares  solutions  become  large.  This  can  be  observed  in  Figure 
3(a)  where  relaxation  of  the  directional  temperatures  is  presented.  When  only  m=2  moments  are 
enforced,  the  use  of  new  velocity-dependent  collision  frequency  results  in  a  reasonable 
approximation  of  the  full  Boltzmann  equation.  However,  as  the  number  of  the  enforced  moments 
increases,  the  reconstruction  of  the  second  moment  becomes  worse. 


(a) 


Time  normalized  to  mean  free  time 


(d) 


Time  normalized  to  mean  free  time 


(g) 


(b) 


(c) 


(e) 


Time  normalized  to  mean  free  time 

(0 


Time  normalized  to  mean  free  time 


(h) 


Time  normalized  to  mean  free  time 
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Time  normalized  to  mean  free  time 


Time  normalized  to  mean  free  time 


Time  normalized  to  mean  free  time 


Time  normalized  to  mean  free  time 


Figure  3.  Relaxation  of  moments  with  kernels  gi,m=(urVi)m,  i  =  1;  2,  m=  2;  3;  4  in  a  mix  of  Maxwellian  streams 

corresponding  to  a  shock  wave  with  Mach  number  6.5. 

Overall,  it  was  concluded  that  the  new  model  is  more  stable  and  allows  the  approximation  of  the 
relaxation  of  the  high  order  moments  with  a  small  number  of  the  basis  functions  in  the 
representation  of  the  collision  frequency.  This  is  an  advantage  in  ID  and  2D  simulations  as 
compared  to  the  old  model.  Indeed,  the  solutions  using  new  model  would  remained  stable,  while 
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the  solution  using  the  old  model  crash  for  the  same  parameters  of  spatial  and  temporal 
discretizations.  Stability  of  the  old  model  can  be  restored,  however,  by  increasing  resolution  in 
physical  space  and  in  time. 

2.4.2  Efficiency  of  the  Evaluation  of  Boltzmann  Collision  Integral  using  Korobov  Nodes 

The  PI  and  the  Graduate  Assistant  evaluated  the  accuracy  and  efficiency  of  the  module  that 
computes  the  Boltzmann  collision  integral  using  the  Korobov  nodes.  It  was  observed  that 
truncation  errors  of  the  Korobov  integration  decrease  faster  than  the  square  root  of  the  number  of 
points  for  smooth  solutions.  It  was  also  observed,  however,  that  accuracy  of  the  schemes  for 
evaluation  of  the  collision  integral  become  acceptable  (less  than  1%)  only  for  large  numbers  of 
Korobov  nodes.  Therefore,  it  was  concluded  that  a  number  of  nodes  comparable  but  smaller  than 
the  number  of  nodes  in  full  Gauss  integration  of  the  collision  integral  is  required  to  achieve  less 
than  1%  l  -error  in  the  collision  integral.  As  a  result,  while  the  method  works  correctly,  it  fails  to 
provide  significant  acceleration  as  compared  to  the  full  deterministic  integration.  The  Pis 
believe,  however,  that  additional  techniques  can  be  devised  that  would  allow  the  use  of  the 
method  in  the  future.  One  technique  that  can  be  tried  is  the  enforcement  of  conservation  laws  in 
the  computed  collision  integral.  The  PI  will  study  this  approach  in  the  future. 

2.4.3  Simulation  of  One  Dimensional  Supersonic  Flow:  the  Normal  Shock  Wave 

The  Pis  applied  the  DGVlib  library  to  the  simulation  of  the  problem  of  a  normal  shock  wave. 
DGVlib  was  merged  with  a  one-dimensional  high-order  Runge-Kutta  discontinuous  Galerkin 
(RKDG)  solver  available  to  the  PI.  The  validity  and  accuracy  of  the  library  was  tested  by 
comparing  to  DSMC  SMILE  solutions  and  experimental  results.  Simulations  were  performed  for 
normal  shock  waves  with  Mach  numbers  1.55,  3.00,  3.8,  6.5,  and  9.0.  The  results  for  shock 
waves  1.55  and  3.8  are  presented  in  Figure  4  and  Figure  5. 

-2 

The  parameters  of  the  downstream  and  upstream  boundary  conditions  are  di=1.0676E-4  kg/m  , 
vi=500.018  m/s,  Ti=300  K  and  d2=1.8990E-4  kg/m3,  v2=28 1.098,  in  the  case  of  the  Mach  1.55 
wave,  and  di=1.0676E-4  kg/m3,  vi=1225.851  m/s,  Ti=300  K  and  d2=3.35386E-4  kg/m3, 
v2=379.  1 322,  in  the  case  of  the  Mach  3.8  wave.  The  gas  is  argon  with  a  specific  heat  ratio  of  5/3. 
Solutions  were  computed  using  four  different  models.  In  the  first  model,  the  ES-BGK  collision 
operator  is  used  with  the  viscosity  law  of  argon  (power  law  exponent  of  .81).  This  solution  is 
marked  “ES-BGK  Argon”  in  Figure  4(a,  b).  In  the  second  model,  the  VD-BGK  collision 
operator  is  used  with  the  relaxation  rate  of  the  directional  temperature  in  the  x-direction 
enforced.  The  relaxation  rates  are  computed  from  the  Boltzmann  collision  operator  with  the 
interval  of  20  mean  free  times  using  a  hard  spheres  collision  model.  The  diameter  of  the  hard 
spheres  is  determined  from  the  variable  soft  spheres  model  (VSS)  so  as  to  match  the  viscosity  of 
the  hard  spheres  gas  to  that  of  viscosity  of  argon.  These  solutions  are  marked  “VD-BGK  Argon” 
in  Figure  4(a,  b).  In  the  third  model,  the  ES-BGK  collision  operator  is  used  with  constants 
corresponding  to  the  viscosity  of  the  hard  spheres  gas  with  the  diameter  of  hard  spheres  set  to 
3.76E-10  m.  These  solutions  are  marked  “ES-BGK”  in  Figure  5(a) — (d).  In  the  fourth  model,  the 
VD-BGK  collision  operator  is  used  in  which  the  relaxation  rate  of  the  directional  temperature  in 
the  x  direction  is  enforced.  The  relaxation  rates  are  determined  by  evaluating  the  Boltzmann 
collision  integral  with  the  interval  of  20  mean  free  times  using  the  hard  spheres  model  with  the 
diameter  of  the  hard  spheres  set  to  3.76E-10  m.  These  solutions  are  marked  “VD-BGK”in  Figure 
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5(a)— (d).  In  Figure  64,  density  profiles  of  the  computed  solutions  are  compared  to  the 
experimental  results  of  Alsmeyer  [7].  In  Figure  5,  solutions  are  compared  to  the  solution  by 
DSMC  solver  SMILE. 


Distance  normalized  to  mean  free  path  Distance  normalized  to  mean  free  path 

Figure  4.  Plots  of  density  and  temperature  profiles  in  solutions  to  the  normal  shock  wave  problem.  Plot  (a) 
corresponds  to  Mach  number  1.55  and  plot  (b)  to  Mach  number  3.8  shock  waves. 

In  the  case  of  the  Mach  1.55  wave,  the  ES-BGK  Argon  solution  is  the  closest  to  the  experimental 
results.  In  the  case  of  the  Mach  3.8  wave,  the  VD-BGK  Argon  is  somewhat  closer  to  the 
experimental  results  than  the  ES-BGK  Argon  solution.  However,  neither  ES-BGK  nor  VD-BGK 
model  achieved  a  good  accuracy.  One  reason  why  the  VD-BGK  model  failed  to  approximate 
experimental  results  may  be  that  the  hard  spheres  collision  model  used  to  determine  the  enforced 
relaxation  rates  may  not  be  the  appropriate  model  to  use  for  simulation  of  Mach  3.8  shock  wave 
in  argon.  Nevertheless,  results  presented  in  Figure  4  illustrate  well  the  fact  that  the  DGVlib 
library  enables  the  user  to  consider  a  whole  range  of  collision  models.  Researchers  will  be  able 
to  select  the  model  that  fits  the  best  in  any  particular  case. 

In  Figure  5  solutions  to  Mach  3.8  shock  wave  obtained  by  ES-BGK  model  and  BGK  model  with 
velocity-dependent  collision  frequency  (VD-BGK)  are  compared  to  the  solution  obtained  by 
DSMC  solver  SMILE.  The  density  profile  of  the  shock  wave  is  presented  in  Figure  5(a).  It  can 
be  seen  that  the  VD-BGK  solution  predicts  density  profile  before  the  shock  better  than  the  ES- 
BGK  solution.  However  both  VD-BGK  and  ES-BGK  miss  the  DSMC  solution  after  the  shock. 
The  velocity  profile  of  the  shock  wave  is  shown  in  Figure  5(b).  It  can  be  seen  that  VD-BGK 
solution  is  significantly  more  accurate  before  the  shock.  After  the  shock  the  ES-BGK  and  VD- 
BGK  solutions  appear  to  be  the  same.  In  Figure  5(c)  the  temperature  profile  of  the  shock  are 
presented  and  in  Figure  5(d)  the  directional  temperature  Tv  is  shown.  Again,  it  can  be  observed 
that  VD-BGK  solutions  are  significantly  more  accurate  in  front  of  the  shock  wave  and  appear 
about  the  same  after  the  shock  wave.  Overall,  it  can  be  concluded  that  the  VD-BGK  model 
yielded  improved  physical  accuracy  of  the  model  kinetic  solutions  as  compared  to  ES-BGK 
model. 
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Figure  5.  Solutions  to  normal  shock  wave  with  Mach  3.8  number. 


2.4.4  Simulation  of  Two-Dimensional  Supersonic  Flow  around  a  Cylinder 

The  main  objective  of  this  task  was  to  implement  the  program  module  for  the  solution  of  the 
right-hand  side  of  the  Boltzmann  and  model  kinetic  equations  developed  at  CSUN  into  a  kinetic 
solver  SMOKE  of  Gimel  Inc,  and  test  the  validity  and  accuracy  of  the  module  through 
comparison  with  benchmark  solutions  obtained  by  the  DSMC  solver  SMILE.  SMOKE  [8]  is  a 
multi-dimensional  finite  volume  solver  to  deterministically  solve  the  BGK  and  ES-BGK  model 
kinetic  equations  on  parallel  computers;  it  is  based  on  conservative  numerical  schemes 
developed  by  L.  Mieussens  [9].  The  kinetic  and  deterministic  nature  of  SMOKE  makes  it  an 
ideal  testing  ground  for  the  developed  program  module.  SMILE  [10]  is  2D/3D  parallel  code  that 
uses  the  DSMC  approach  for  the  statistical  solution  of  the  full  Boltzmann  equation;  it  has  been 
extensively  validated  for  various  problems  through  comparison  with  experimental  data  and  other 
numerical  solutions. 


2.4.4.1  Flow  Conditions 

The  problem  under  consideration  is  supersonic  flow  over  a  two-dimensional  cylinder.  The  free 
stream  gas  is  nitrogen,  mass  density  is  le-6  kg/m3,  temperature  250  K,  and  velocity  1,000  m/s. 
This  corresponds  to  a  ~M=3,  Kn=0.1  flow  based  on  a  0.5  m  cylinder  diameter.  Exponential 
viscosity  -  temperature  dependence  is  used,  with  the  exponent  values  of  0.5  (hard  sphere  model) 
and  0.75  used  in  the  computations.  The  fully  diffuse  reflection  on  cylinder  wall  is  used,  with  the 
wall  temperature  set  first  to  300  K  and  then  to  500  K. 
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2.4.4.2  Grid  Convergence  Study 

Prior  to  the  validation  study  and  accuracy  analysis,  the  spatial  and  velocity  grid  resolutions 
necessary  for  obtaining  grid  independent  results  were  determined.  The  original  SMOKE  was 
used  for  this  purpose.  The  number  of  spatial  cells  varied  from  5,000  to  over  25,000,  and  the 
number  of  velocity  bins  in  each  of  the  three  directions  varied  from  20  to  60.  An  example  of  the 
convergence  results  is  shown  in  Figure  6,  where  the  gas  macroparameters  obtained  for  two 
different  velocity  grid  resolutions  are  presented. 


The  example  shows  that  even  a  20x20x20  velocity  grid  resolution  provides  high  accuracy  in 
front  of  the  cylinder,  although  there  is  some  small  difference  in  the  wake  flow  where  there  are 
clearly  more  than  20  velocity  points  needed.  The  results  showed  that  a  velocity  grid  of 
40x40x40,  along  with  5,000  spatial  cells,  provide  adequate  resolution  as  further  refinement  does 
not  change  the  macroparameters. 
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Figure  6.  Density  (left)  and  temperature  (right)  fields  around  a  cylinder  obtained  for  two  different  velocity  grids 
(20x20x20,  top,  and  40x40x40,  bottom,  for  X-Y-Z  directions).  TwaM=300  K. 


2. 4.4. 3  Adding  ES-BGK  Capability 

The  incorporation  of  the  developed  collision  module  into  SMOKE  allowed  a  direct  comparison 
of  the  original  and  extended  SMOKE  (called  SMOKE-DGV  hereafter)  in  terms  of  the  solution  of 
the  ES-BGK  equation.  Comparison  of  the  two  codes,  SMOKE  and  SMOKE-DGV,  is  presented 
in  Figure  7.  One  can  see  that  there  is  quite  reasonable  agreement  between  the  implementations, 
even  though  different  methods  of  the  solution  of  the  ES-BGK  equation  are  used.  Good 
agreement  between  the  codes  provides  verification  of  the  implementation  of  the  DGV  extension. 

As  expected,  the  solution  of  the  ES-BGK  equation  is  noticeably  different  from  the  Boltzmann 
equation;  consistent  with  similar  comparisons  published  in  the  literature  rill,  the  ES-BGK 
solution  is  characterized  by  a  somewhat  wider  shock  front.  The  difference  in  the  wake  region  is 
less  significant  than  that  in  the  shock  since  the  non-equilibrium  is  less  pronounced  in  that  region. 
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Note  that  an  ES-BGK  solution  is  expected  to  be  more  accurate  for  Maxwellian  molecules,  and 
less  accurate  for  hard  spheres.  This  is  illustrated  in  Figure  8,  where  the  temperature  fields 
obtained  by  SMILE  and  SMOKE-DGV  are  shown.  It  can  be  seen  that,  although  there  is  good 
agreement  in  the  wake  region,  the  bow  shock  is  much  thicker  in  the  ES-BGK  solution. 
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Figure  7.  Velocity  in  X  direction  obtained  by  the  solution  of  the  ES-BGK  equation  with  the  original  and  extended 

SMOKE. 
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Figure  8.  Comparison  of  the  Boltzmann  and  ES-BGK  solutions  for  Maxwellian  molecules  (left)  and  hard  sphere 

molecules  (right),  TwaM=500  K. 
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2.4.4.4  BGK  with  Velocity-Dependent  Collision  Frequency:  Correction  through  Boltzmann 
Equation 

Let  us  now  examine  the  effect  of  using  a  velocity-dependent  collision  frequency  for  the  ES-BGK 
equation,  obtained  through  the  evaluation  of  the  full  Boltzmann  collision  integral.  The  gas 
properties  obtained  for  the  hard  sphere  model  are  presented  in  Figure  9.  Here,  the  bottom  halves 
show  the  results  computed  with  extended  SMOKE  ES-BGK  models  (SMOKE-DGV)  and  the 
model  with  velocity-dependent  collision  frequency  (SMOKE- vDGV).  It  can  be  seen  that  the 
correction  of  the  collision  frequency  allows  for  much  closer  agreement  with  the  solution  of  the 
Boltzmann  equation.  Note  that  the  correction  was  conducted  about  every  500-th  time  step,  and 
the  total  computational  time  increased  by  approximately  50%.  Still,  further  refinement  through 
the  increase  in  the  resolution  of  the  velocity  grid  during  the  evaluation  of  the  Boltzmann 
collision  integral  (currently  at  15x15x15)  requires  significantly  larger  computer  resources,  and 
thus  efficient  parallelization. 
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Figure  9.  Velocity  field  (left)  and  temperature  (right)  obtained  by  the  solution  of  the  Boltzmann  equation  (top 

halves)  and  velocity-dependent  ES-BGK  (bottom  halves). 

Also  note  that  a  non-equilibrium  dependent  switch  has  been  used  between  the  regular  and 
velocity-dependent  modes  of  the  ES-BGK  equation.  In  regions  where  non-equilibrium  is 
significant,  the  velocity-dependent  mode  is  used,  while  in  the  rest  of  the  domain,  the 
conventional  ES-BGK  equation  is  solved.  The  degree  of  non-equilibrium  is  evaluated  by  the 
A/=!/-/e9l// metric,  where/ and  feq  are  the  local  computed  and  corresponding  equilibrium 
distribution  functions.  The  values  of  this  metric  are  given  in  Figure  10.  The  velocity-dependent 
collision  frequency  was  used  in  regions  where  A/  >0.05. 
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Figure  10.  Non-equilibrium  metric  used  for  switching  between  conventional  and  velocity-dependent  frequency 

ES-BGK  equations. 


3.  SUMMARY 

New  modules  for  evaluating  the  kinetic  collision  operator  using  BGK  model  with  velocity- 
dependent  collision  frequency  (VD-BGK)  and  for  evaluating  the  Boltzmann  collision  integral 
using  quasi-stochastic  integration  have  been  implemented.  Benchmark  solutions  have  been 
computed  to  test  the  VD-BGK  model.  Solutions  were  compared  to  DSMC  solutions  and 
experimental  results. 

In  zero,  one-  and  two-dimensional  simulations  it  was  observed  that  the  Boltzmann  equation 
based  correction  to  the  collision  frequency  improves  agreement  with  the  solution  of  the 
Boltzmann  equation  as  compared  to  the  classical  ES-BGK  model  equations.  One  difficulty  that 
persists  is  the  fact  that  the  VD-BGK  model  makes  the  kinetic  equation  considerably  stiffer  when 
high  order  moments  are  enforced.  The  new  least-squares  approach  helps  to  reduce  the  stiffness 
significantly,  however,  different  forms  of  velocity-dependent  collision  frequency  may  need  to  be 
explored  to  remedy  the  excess  stiffness  entirely.  The  need  to  evaluate  the  Boltzmann  collision 
operator  to  estimate  relaxation  rates  of  the  moments  in  the  VD-BGK  model  is  limiting  efficiency 
of  the  model,  especially  for  two  dimensional  problems.  The  introduction  of  coarser  secondary 
grid  in  DGVlib  for  such  estimates  did  improve  the  efficiency.  However,  shock  waves  with  high 
Mach  numbers  still  require  large  secondary  grids.  A  possible  remedy  to  that  may  be  a  direct 
computation  of  the  moments  of  the  collision  operator  bypassing  the  evaluation  of  the  collision 
operator  itself.  These  evaluations  are  orders  of  magnitude  faster  than  evaluations  of  the  full 
collision  integral,  but  require  significantly  more  memory  for  storing  pre-computed  moment 
kernels.  Formulating  and  implementing  this  approach  will  be  Pis  future  work. 

Korobov  integration  of  the  collision  integral  proved  to  be  more  evasive  than  originally  thought. 

It  was  observed  that  Korobov  quadratures  do  converge  to  the  desired  integral.  However, 
practically  acceptable  levels  of  errors  are  achieved  only  for  numbers  of  integration  points 
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comparable  to  those  in  full  deterministic  evaluation  of  the  collision  integral.  The  Pis  hope, 
however,  that  additional  remedies  can  be  introduced  to  the  method,  e.g.,  enforcement  of 
conservation  laws  in  the  computed  collision  integral.  These  remedies  would  allow  use  of  smaller 
number  of  quadrature  points.  The  proposed  future  modification  of  the  methods  could  be 
implemented  in  DGVlib  with  a  minor  programming  effort.  Therefore,  DGVlib  represents  a 
flexible  software  tool  that  can  be  easily  modified  to  improve  the  models. 
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